/*
  tlc_s05: template light curve (TLC) fitting program to fit sparse JHK-band Cepheid light curves with TLC 
           from Soszynski et al (2005, PASP 117:823, hereafter s05)

  Input: two required inputs
    -infile <input_filename> -- input ASCII datafile with either of the following two formats
                                        time  mag  magerr
                                           --  OR --
                                       phase  mag  magerr
    -band <J,H,K>            -- set the fitting of template light curves for either J, H or K band

  Example Usage:
    ./tlc_s05 -infile AAGem_J.dat -band J
    
  Optional Inputs:
     -period <PERIOD> -- use the input period to phase the data if input file is in [time mag magerr] format
     -t0 <TIME_0>     -- for phasing the input file with this epoch
     -v               -- use the V-band based template light curves (see Soszynski et al for more detail)
     -i               -- use the I-band based template light curves (see Soszynski et al for more detail)
   [Note: the default is using the average of the V- and I-band based template light curves]

  Output on Screen: 
    input_filename      intensity_mean     fitted_amplitude     [chi2]

  Output ASCII Files: 
     [input_filename].phdata -- phased input data file in the form of: phase mag mag_err
     [input_filename].lcurve -- fitted template light curve in the form of: phase mag

  Required external package: GSL (GNU Scientific Library, http://www.gnu.org/software/gsl/)
    Note: the GSL is assumed to be installed in standard directory (e.g., /usr/local/lib). 
          If not link to the GSL library needs to be changed when compiling the code

  Compile: to compile the code, type the following 
    gcc -o tlc_s05 tlc_s05.c -L/usr/local/lib -lgsl -lgslcblas -lm

  Reference:
    Some codes, especially the GSL part, has been referred from http://lists.gnu.org/archive/html/help-gsl/2009-12/msg00034.html

  Author: Chow-Choong Ngeow (National Central University)
  
  Date first created: 12 Sept 2014

  Citation: Please cite the following if you use this program
    Ngeow & Kanbur, 2014, Proceeding of APRIM2014 [TBD]

  Copyright 2014 Chow-Choong Ngeow

  License: GNU GENERAL PUBLIC LICENSE Version 3 

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

  Revision history:
    Version 1.0 -- initial version

*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

/* --- headers for GSL --- */
#include "gsl/gsl_multimin.h"
#include "gsl/gsl_statistics.h"
#include "gsl/gsl_sort.h"
#include "gsl/gsl_sort_double.h"
#include "gsl/gsl_sort_vector.h"

/* --- various definition --- */
#define VERSION "1.0"

#define TEMPLATE_AVE 0
#define TEMPLATE_V   1
#define TEMPLATE_I   2

#define BAND_J 1
#define BAND_H 2
#define BAND_K 3

#define N_ORDER 7

/* --- data structure --- */
struct data {
  int         n;   
  gsl_vector *x;   
  gsl_vector *y;   
  gsl_vector *w;   
};

/* --- global variables --- */
double amp_template[N_ORDER+1],phi_template[N_ORDER+1];

/* --- main subroutines header --- */
void get_template(int TEMPLATE, int BAND);
double phase_adjust(double inphase);
double recalculate_chi2(int N, char *filename);
double fit_templateLC(int N, double phase[], double mag[], double err[], double werr[], double mean_in, double ampl_in, double delta_phase, double *mean_fit, double *ampl_fit);

void print_usage(char *code)
{
  printf("Usage: %s -infile <Input.dat> -band <J, H or K>\n",code);
  printf("\t Option:\n");
  printf("\t\t -i                Use I-band based template light curves\n");
  printf("\t\t -v                Use V-band based template light curves\n");
  printf("\t\t -p    <#>         Period in days\n");
  printf("\t\t -t0   <#>         Value of the epoch (default: 0.0)\n");
  printf("\t\t -version          Print current version\n");
  
  printf("\n");
  printf("  Input data file format: [time mag err]  OR  [phase mag err]\n");
  printf("\n");
}

double my_f (const gsl_vector *v, void *params) 
{ 
  struct data *p = (struct data *) params;

  int i,j;
  double chi2;
  double xi,yi,wi,y_calc;
  double m = gsl_vector_get (v, 0);
  double amp = gsl_vector_get (v, 1);
  double tsum;
  double *ys;

  ys = (double *)calloc(p->n,sizeof(double ));

  for(i=0;i<p->n;i++) {

    xi = gsl_vector_get(p->x, i);

    tsum=0.0;
    for(j=1;j<=N_ORDER;j++)
      tsum += amp_template[j] * cos(2.0*M_PI*j*xi + phi_template[j]);
    
    ys[i] = tsum;
  }

  chi2=0.0;
  for (i=0;i<p->n;i++) {

    yi = gsl_vector_get(p->y, i);    
    wi = gsl_vector_get(p->w, i);
    
    y_calc = m + fabs(amp)*ys[i];
    chi2 += wi * ((yi - y_calc) * (yi - y_calc));
  }
  
  free(ys);

  return chi2;
}

int main(int argc, char *argv[])
{
  int i,j,N;
  int input_band=0,flag_in=0,flag_p=0,flag_template=TEMPLATE_AVE;
  char infile[3000],outphfile[3000],outlcfile[3000],line[2000];
  double period,t0=0.0;
  double t_time,t_mag,t_err,t_phase;
  double *phase,*mag,*err,*werr;
  double chi2,chi2_min,delta_phase,dphase_best;
  double mean_fit,mean_best,mean_in;
  double ampl_fit,ampl_best,ampl_in;
  double sum,intmean,UI;
  double *out_phase,*out_mag;

  int get_input();				
  double get_input_dbl();
  void print_usage();
  void get_template();

  FILE *fin,*fout;

  /* --- print out usage --- */
  if(argc<2) {
    print_usage(argv[0]);
    exit(0);
  }  

  /* --- input command line argument --- */
  for(i=1;i<argc;i++) {
    
    if(!strcmp("-infile",argv[i]))
      flag_in=get_input("infile",argv[i+1],infile);
    
    if(!strcmp("-band",argv[i])) {
      if(!strcmp("J",argv[i+1]) || !strcmp("j",argv[i+1])) input_band = BAND_J;
      else if(!strcmp("H",argv[i+1]) || !strcmp("h",argv[i+1])) input_band = BAND_H;
      else if(!strcmp("K",argv[i+1]) || !strcmp("k",argv[i+1])) input_band = BAND_K;
      else {
	printf(" Error: Input for option -band is not J, H or K, abort\n");
	exit(0);
      }
    }
    
    /* --- options --- */
    if(!strcmp("-v",argv[i])) flag_template=TEMPLATE_V;

    if(!strcmp("-i",argv[i])) flag_template=TEMPLATE_I;
    
    if(!strcmp("-period",argv[i])) {
      period=get_input_dbl("period",argv[i+1]);
      flag_p=1;
    }
    
    if(!strcmp("-t0",argv[i])) 
      t0=get_input_dbl("t0",argv[i+1]);
    
    
    if(!strcmp("-version",argv[i])) { 
      printf("%s: Version %s\n",argv[0], VERSION);
      exit(0);
    }

    /* --- print out helps --- */
    if(!strcmp("-h",argv[i]) || !strcmp("--h",argv[i]) || !strcmp("-help",argv[i]) || !strcmp("--help",argv[i])) {
      print_usage(argv[0]);
      exit(0);
    }  
  }
  
  /* --- check input file argument is given --- */
  if(!flag_in) {
    printf(" Error: required -infile <data.file> not given \n");
    exit(0);
  }

  /* --- setup the outputs files --- */
  sprintf(outphfile,"%s.phdata",infile);
  sprintf(outlcfile,"%s.lcurve",infile);

  /* --- count the lines in input file --- */
  N=0;
  fin=fopen(infile,"r"); 
  if(fin==NULL) {
    printf(" Error: cannot open file %s\n",infile);
    exit(0);
  }
  else {
    while(fgets(line,2000,fin)!=NULL)
      N++;
    fclose(fin);
  }

  /* --- can only use for N>=3 --- */
  if(N<3) {
    printf(" Warning: %s has less than 3 data points, abort\n",infile);
    exit(0);
  }
  
  /* --- memory allocation, input data and output the phased light curve --- */
  phase=(double *)calloc(N,sizeof(double ));
  mag=(double *)calloc(N,sizeof(double ));
  err=(double *)calloc(N,sizeof(double ));
  werr=(double *)calloc(N,sizeof(double ));

  fin=fopen(infile,"r"); fout=fopen(outphfile,"w");
  for(i=0;i<N;i++) {

    fscanf(fin,"%lg %lg %lg",&t_time,&t_mag,&t_err);

    mag[i]=t_mag;
    err[i]=t_err;
    werr[i]=1.0/(t_err*t_err);

    /* --- convert to phase if period is given --- */
    if(flag_p==1)
      phase[i]=(t_time-t0)/period - (int)((t_time-t0)/period);
    else
      phase[i]=t_time;

    /* --- make sure the phase is positive --- */
    if(phase[i]<0.0)
      phase[i]+=1.0;
    
    /* --- output the phased light curves --- */
    fprintf(fout,"%2.3f\t%2.3f\t%2.3f\n",phase[i],mag[i],err[i]);
  }
  fclose(fin); fclose(fout);

  /* --- setup the template light curves --- */
  get_template(flag_template, input_band);
  
  /* ======================================== */
  /*                                          */
  /*    TEMPLATE LIGHT CURVE FITTING HERE     */
  /*                                          */
  /* ======================================== */

  /* --- setup the input mean mag and amplitude --- */
  mean_in=gsl_stats_wmean(werr,1,mag,1,N);
  ampl_in=0.4;

  /* --- loop over delta_phase to fit the template light curve --- */
  chi2_min=90000000.0;
  for(i=0;i<=1000;i++) {
     
    delta_phase=0.001*i;

    chi2=fit_templateLC(N, phase, mag, err, werr, mean_in, ampl_in, delta_phase, &mean_fit, &ampl_fit);

    if(chi2<chi2_min) {
      chi2_min=chi2;
      mean_best=mean_fit;
      ampl_best=ampl_fit;
      dphase_best=delta_phase;
    }
  }

  ampl_best = fabs(ampl_best);

  /* --- output best fit phased template light curve --- */
  out_phase=(double *)calloc(1000,sizeof(double ));
  out_mag=(double *)calloc(1000,sizeof(double ));
  
  for(i=0;i<1000;i++) {
    
    t_phase = 0.001*i;
    
    t_mag=mean_best;
    for(j=1;j<=N_ORDER;j++)
      t_mag += ampl_best*(amp_template[j] * cos(2.0*M_PI*j*t_phase + phi_template[j]));
    
    out_phase[i]=phase_adjust(t_phase+(1.0-dphase_best));
    out_mag[i]=t_mag;
  }
  
  gsl_sort2(out_phase, 1, out_mag, 1, 1000);

  fout=fopen(outlcfile,"w");
  for(i=0;i<1000;i++)     
    fprintf(fout,"%2.3f\t%2.3f\n",out_phase[i],out_mag[i]);
  fclose(fout);

  /* --- calculate the intensity mean from best-fit template light curve --- */
  sum=0.0;
  fin=fopen(outlcfile,"r");
  for(i=0;i<1000;i++) {
    fscanf(fin,"%lg %lg",&t_phase,&t_mag);
    sum+= pow(10.0,-0.4*t_mag);
  }
  fclose(fin);

  intmean = -2.5*log10(sum/1000.0);

  /* --- recalculate the chi2 value --- */
  chi2_min = recalculate_chi2(N,infile);

  /* --- print out results --- */
  printf("%s\t%2.3f\t%2.3f \t[chi2= %2.4f]\n",infile,intmean,ampl_best,chi2_min);

  /* --- free memory --- */
  free(phase); free(mag); free(err); free(werr);
  free(out_phase); free(out_mag);

  return 0;
}

/* --- subroutine to setup template light curve --- */
void get_template(int TEMPLATE, int BAND)
{
  // J band
  if(BAND == 1)
    {
      switch(TEMPLATE) {
	// V-band based
      case 1:
	amp_template[1] = 0.432; phi_template[1] = 1.734;
	amp_template[2] = 0.109; phi_template[2] = 2.747;
	amp_template[3] = 0.058; phi_template[3] = 3.079;
	amp_template[4] = 0.039; phi_template[4] = 3.671;
	amp_template[5] = 0.025; phi_template[5] = 4.253;
	amp_template[6] = 0.016; phi_template[6] = 4.460;
	amp_template[7] = 0.010; phi_template[7] = 5.429;
	break;
      
	// I-band based
      case 2:
	amp_template[1] = 0.432; phi_template[1] = 1.863;
	amp_template[2] = 0.110; phi_template[2] = 2.901;
	amp_template[3] = 0.060; phi_template[3] = 3.324;
	amp_template[4] = 0.038; phi_template[4] = 3.934;
	amp_template[5] = 0.029; phi_template[5] = 4.678;
	amp_template[6] = 0.015; phi_template[6] = 5.164;
	amp_template[7] = 0.010; phi_template[7] = 5.851;
	break;
      
	// average of VI-based
      default:
	amp_template[1] = 0.432; phi_template[1] = 1.799;
	amp_template[2] = 0.110; phi_template[2] = 2.824;
	amp_template[3] = 0.059; phi_template[3] = 3.202;
	amp_template[4] = 0.039; phi_template[4] = 3.803;
	amp_template[5] = 0.027; phi_template[5] = 4.466;
	amp_template[6] = 0.016; phi_template[6] = 4.812;
	amp_template[7] = 0.010; phi_template[7] = 5.640;
      }
    }
  
  // H band
  if(BAND == 2)
    {
      switch(TEMPLATE) {
	// V-band based
      case 1:
	amp_template[1] = 0.438; phi_template[1] = 1.262;
	amp_template[2] = 0.089; phi_template[2] = 2.303;
	amp_template[3] = 0.039; phi_template[3] = 2.913;
	amp_template[4] = 0.019; phi_template[4] = 3.394;
	amp_template[5] = 0.019; phi_template[5] = 3.318;
	amp_template[6] = 0.018; phi_template[6] = 3.972;
	amp_template[7] = 0.011; phi_template[7] = 4.727;
	break;
      
	// I-band based
      case 2:
	amp_template[1] = 0.433; phi_template[1] = 1.397;
	amp_template[2] = 0.094; phi_template[2] = 2.543;
	amp_template[3] = 0.046; phi_template[3] = 3.266;
	amp_template[4] = 0.023; phi_template[4] = 3.886;
	amp_template[5] = 0.022; phi_template[5] = 3.813;
	amp_template[6] = 0.021; phi_template[6] = 4.380;
	amp_template[7] = 0.016; phi_template[7] = 5.245;
	break;
      
	// average of VI-based
      default:
	amp_template[1] = 0.436; phi_template[1] = 1.330;
	amp_template[2] = 0.092; phi_template[2] = 2.423;
	amp_template[3] = 0.043; phi_template[3] = 3.090;
	amp_template[4] = 0.021; phi_template[4] = 3.630;
	amp_template[5] = 0.021; phi_template[5] = 3.566;
	amp_template[6] = 0.020; phi_template[6] = 4.176;
	amp_template[7] = 0.014; phi_template[7] = 4.986;
      }
    }

  // K band
  if(BAND == 3)
    {
      switch(TEMPLATE) {
	// V-band based
      case 1:
	amp_template[1] = 0.440; phi_template[1] = 1.201;
	amp_template[2] = 0.082; phi_template[2] = 2.174;
	amp_template[3] = 0.035; phi_template[3] = 2.757;
	amp_template[4] = 0.021; phi_template[4] = 3.448;
	amp_template[5] = 0.012; phi_template[5] = 3.406;
	amp_template[6] = 0.013; phi_template[6] = 3.658;
	amp_template[7] = 0.010; phi_template[7] = 4.719;
	break;
	
	// I-band based
      case 2:
	amp_template[1] = 0.434; phi_template[1] = 1.333;
	amp_template[2] = 0.085; phi_template[2] = 2.419;
	amp_template[3] = 0.043; phi_template[3] = 3.097;
	amp_template[4] = 0.020; phi_template[4] = 3.794;
	amp_template[5] = 0.016; phi_template[5] = 3.765;
	amp_template[6] = 0.015; phi_template[6] = 4.475;
	amp_template[7] = 0.013; phi_template[7] = 5.149;
	break;
      
	// average of VI-based
      default:
	amp_template[1] = 0.437; phi_template[1] = 1.267;
	amp_template[2] = 0.084; phi_template[2] = 2.297;
	amp_template[3] = 0.039; phi_template[3] = 2.927;
	amp_template[4] = 0.021; phi_template[4] = 3.621;
	amp_template[5] = 0.014; phi_template[5] = 3.586;
	amp_template[6] = 0.014; phi_template[6] = 4.067;
	amp_template[7] = 0.012; phi_template[7] = 4.934;
      }     
    }
}


/* --- fitting template light curve and return chi2 value --- */
double fit_templateLC(int N, double phase[], double mag[], double err[], double werr[], double mean_in, double ampl_in, double delta_phase, double *mean_fit, double *ampl_fit)
{
  int i,j,nparas=2;
  double phtem,magtem;
  double chi2fit;
  double *x,*y,*w;
  struct data data_input;

  /* --- memory allocation for x,y,w vectors and setup data arrays for GSL --- */ 
  x=(double *)calloc(N,sizeof(double ));
  y=(double *)calloc(N,sizeof(double ));
  w=(double *)calloc(N,sizeof(double ));

  data_input.n = N;
  data_input.x = gsl_vector_alloc(data_input.n); 
  data_input.y = gsl_vector_alloc(data_input.n); 
  data_input.w = gsl_vector_alloc(data_input.n); 
  for (i=0;i<data_input.n;i++) { 
    gsl_vector_set (data_input.x, i, phase_adjust(phase[i]+delta_phase));
    gsl_vector_set (data_input.y, i, mag[i]);
    gsl_vector_set (data_input.w, i, werr[i]);
  }

  /* --- initial guess --- */
  gsl_vector *paras;
  paras = gsl_vector_alloc (nparas);
  gsl_vector_set (paras, 0, mean_in);
  gsl_vector_set (paras, 1, ampl_in);

  /* --- initialise algorithm parameters --- */
  size_t iter = 0;
  int status;
  double size;

  /* --- set initial STEP SIZES to 0.001 --- */
  gsl_vector *ss;
  ss = gsl_vector_alloc (nparas); 
  gsl_vector_set_all (ss, 0.001);
  
  const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
  gsl_multimin_fminimizer *s = NULL;
  gsl_multimin_function minex_func;
  

  /* --- initialize method and iterate --- */
  minex_func.n = nparas;
  minex_func.f = my_f;
  minex_func.params = &data_input;
  
  s = gsl_multimin_fminimizer_alloc (T, nparas);
  gsl_multimin_fminimizer_set (s, &minex_func, paras, ss); 
  
  do {
    iter++;
    status = gsl_multimin_fminimizer_iterate(s);
    
    if (status)
      break;
    
    size = gsl_multimin_fminimizer_size (s);
    status = gsl_multimin_test_size (size, 1e-5);
  } 
  while (status == GSL_CONTINUE && iter < 1000);

  if (status == GSL_SUCCESS) {

    *mean_fit = gsl_vector_get (s->x, 0);
    *ampl_fit = gsl_vector_get (s->x, 1);
    chi2fit = s->fval;
  }

  /* --- free memory and return chi2 --- */
  gsl_multimin_fminimizer_free (s);  
  gsl_vector_free(data_input.x);
  gsl_vector_free(data_input.y);
  gsl_vector_free(data_input.w);
  gsl_vector_free(paras);
  gsl_vector_free(ss);
  
  free(x); free(y); free(w);

  return chi2fit;
}

/* ---recalculate the chi2 based on two output files --- */
double recalculate_chi2(int N, char *filename)
{
  char infile[1000];
  int i,j;
  double *phs_tt,*mag_tt,*err_tt;
  double *phs_tl,*mag_tl;
  double sum,chi2;

  FILE *fin;

  /* --- memory allocation --- */
  phs_tt=(double *)calloc(N,sizeof(double ));
  mag_tt=(double *)calloc(N,sizeof(double ));
  err_tt=(double *)calloc(N,sizeof(double ));
  phs_tl=(double *)calloc(1000,sizeof(double ));
  mag_tl=(double *)calloc(1000,sizeof(double ));

  /* --- input data -- */
  sprintf(infile,"%s.phdata",filename);
  fin=fopen(infile,"r");
  for(i=0;i<N;i++)
    fscanf(fin,"%lg %lg %lg",&phs_tt[i],&mag_tt[i],&err_tt[i]);
  fclose(fin);

  sprintf(infile,"%s.lcurve",filename);
  fin=fopen(infile,"r");
  for(i=0;i<1000;i++)
    fscanf(fin,"%lg %lg",&phs_tl[i],&mag_tl[i]);
  fclose(fin);

  /* --- calculate chi2 --- */
  sum = 0.0;
  for(i=0;i<1000;i++) {

    for(j=0;j<N;j++) {

      if(phs_tt[j] == phs_tl[i]) 
	sum += (mag_tt[j] - mag_tl[i]) * (mag_tt[j] - mag_tl[i]) / (err_tt[j]*err_tt[j]); 
    }
  }

  chi2 = sum;
	  
  /* --- free memory and return chi2 --- */
  free(phs_tt); free(mag_tt); free(err_tt);
  free(phs_tl); free(mag_tl);

  return chi2;
}

/* --- adjust the phase value --- */
double phase_adjust(double inphase)
{
  double newphase;

  if(inphase>=1.0) newphase=inphase-1.0;
  else if(inphase<=0.0) newphase=inphase+1.0;
  else newphase=inphase;

  return newphase;
}

/* --- various subroutines for command-line inputs --- */
int get_input(char *option, char *input, char *output)
{
  if(input!=NULL) {
    sprintf(output,"%s",input);
    return (1);
  }
  else {
    printf(" Error: Input for option -%s not given, abort\n",option);
    exit(0);
  }
}

double get_input_dbl(char *option, char *input)
{
  double var;
  if(input!=NULL) 
    sscanf(input,"%lg",&var);
  else {
    printf(" Error: Input for option -%s not given, abort\n",option);
    exit(0);
  }
  return (var);
}


#undef VERSION

#undef TEMPLATE_AVE 
#undef TEMPLATE_V   
#undef TEMPLATE_I   

#undef BAND_J
#undef BAND_H 
#undef BAND_K 

#undef N_ORDER
